function [S0,P0,P1] = kf_SWR(YS1,YS2,Q,R,Sm,rho,SI,PI,TT);


% This file performs the forward kalman filter recursions for the Stock-Watson-Romer model.  
% YS1,YS2 are the noisy and clean measures of inflation, respectively
% Q,R are the SW state innovation variances
% Sm is the standard deviation of the measurement error
% rho is the ar parameter of the measurement error
% SI,PI are the initial values for the recursions, S(1|0) and P(1|0)
% TT = [T1,T2,T] 

T1 = TT(1); %1947
T2 = TT(2); %1990
T = TT(3);  %2013

% States are ordered \pi, \mu, m
S0 = zeros(3,T); % current estimate of the state, S(t|t)
S1 = zeros(3,T); % one-step ahead estimate of the state, S(t+1|t)
P0 = zeros(3,3,T); % current estimate of the covariance matrix, P(t|t)
P1 = zeros(3,3,T); % one-step ahead covariance matrix, P(t+1|t)

% constant parameters
A = [0 1 0; 0 1 0; 0 0 rho];

% date 1
C = [1 0 1];
D = 0;
y10 = C*SI; % E(y(t|t-1)
V10 = C*PI*C' + D*D'; % V(y(t|t-1)
S0(:,1) = SI + PI*(C'/V10)*( YS1(:,1) - y10 ); % E(S(t|t))
P0(:,:,1) = PI - PI*(C'/V10)*C*PI; % V(S(t|t))
S1(:,1) = A*S0(:,1); % E(S(t+1|t)
B = [R(2)^.5 Q(2)^.5 0; 0 Q(2)^.5 0; 0 0 Sm];
P1(:,:,1) = A*P0(:,:,1)*A' + B*B'; % V(S(t+1|t)

% Iterating through 1947 (single noisy observation)
for i = 2:T1,
  y10 = C*S1(:,i-1); % E(y(t|t-1)
  V10 = C*P1(:,:,i-1)*C' + D*D'; % V(y(t|t-1)
  S0(:,i) = S1(:,i-1) + P1(:,:,i-1)*(C'/V10)*( YS1(:,i) - y10 ); % E(S(t|t))
  P0(:,:,i) = P1(:,:,i-1) - P1(:,:,i-1)*(C'/V10)*C*P1(:,:,i-1); % V(S(t|t))
  S1(:,i) = A*S0(:,i); % E(S(t+1|t)
  B = [R(i+1)^.5 Q(i+1)^.5 0; 0 Q(i+1)^.5 0; 0 0 Sm];
  P1(:,:,i) = A*P0(:,:,i)*A' + B*B'; % V(S(t+1|t)
end

% Iterating 1948-1990 (two observations)
C = [1 0 1; 1 0 0];
D = sqrt(eps)*eye(2); % D set to small values to simplify programming
for i = T1+1:T2,
  y10 = C*S1(:,i-1); % E(y(t|t-1)
  V10 = C*P1(:,:,i-1)*C' + D*D'; % V(y(t|t-1)
  S0(:,i) = S1(:,i-1) + P1(:,:,i-1)*(C'/V10)*( [YS1(:,i); YS2(:,i)] - y10 ); % E(S(t|t))
  P0(:,:,i) = P1(:,:,i-1) - P1(:,:,i-1)*(C'/V10)*C*P1(:,:,i-1); % V(S(t|t))
  S1(:,i) = A*S0(:,i); % E(S(t+1|t)
  B = [R(i+1)^.5 Q(i+1)^.5 0; 0 Q(i+1)^.5 0; 0 0 Sm];
  P1(:,:,i) = A*P0(:,:,i)*A' + B*B'; % V(S(t+1|t)
end

% Iterating 1991-2013 (single clean observation)
C = [1 0 0];
D = sqrt(eps);
for i = T2+1:T,
  y10 = C*S1(:,i-1); % E(y(t|t-1)
  V10 = C*P1(:,:,i-1)*C' + D*D'; % V(y(t|t-1)
  S0(:,i) = S1(:,i-1) + P1(:,:,i-1)*(C'/V10)*( YS2(:,i) - y10 ); % E(S(t|t))
  P0(:,:,i) = P1(:,:,i-1) - P1(:,:,i-1)*(C'/V10)*C*P1(:,:,i-1); % V(S(t|t))
  S1(:,i) = A*S0(:,i); % E(S(t+1|t)
  B = [R(i+1)^.5 Q(i+1)^.5 0; 0 Q(i+1)^.5 0; 0 0 Sm];
  P1(:,:,i) = A*P0(:,:,i)*A' + B*B'; % V(S(t+1|t)
end
